%% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function varargout = ModelPoles(cfg,exvars,depvars,indepvars,r,time,varargin)
%% Destinguish between calibration or validation and specific type of run
if ~strcmp(cfg.ModeSet,'modes4')
    error ('The Poles model requires modes4 for the modeset')
end

for t=time
    for n = 2:2:length(varargin)
        eval( [varargin{n-1} '(r,1,:)=varargin{n};']);
    end
    % Get rid of walking and biking modes
    modes  =  indepvars.UseModes(r,1,:)& ~(isequal(cfg.modes,'UM')|isequal(cfg.modes,'Bike')|isequal(cfg.modes,'Walk') );
    exvars.TravePrice = repmat(exvars.GasolinePrice,[1 1 4]);
    
if t == 1 || t == 2
    %% POLES model has to initialize values before running
    depvars.TravelDemand(r,t,modes)  = exvars.TravelDemand (r,t,modes);
    depvars.EquipmentRate(r,t,3)     = exvars.EquipmentRate(r,t);
    depvars.VTravelDemand(r,t,3)     = depvars.TravelDemand(r,t,3)/depvars.EquipmentRate(r,t,3);
    
else
    %% Remaining runs
    % Cars
    depvars.IncomeElasticityER(r,t,3) = indepvars.IncomeElasticityER(r,1,1) .* exp((depvars.EquipmentRate(r,t-1,3)/...
    (depvars.EquipmentRate(r,t-1,3) - exvars.ERSaturation(r,1,1)))*0.05);

    depvars.IncomeElasticity(r,t,4) = indepvars.IncomeElasticity(r,1,4) .* exp((depvars.TravelDemand(r,t-1,4)./...
    (depvars.TravelDemand(r,t-1,4) -  indepvars.Saturation(r,1,4))).*0.05);
    
    depvars.EquipmentRate(r,t,3) = depvars.EquipmentRate(r,t-1,3) *...
        max([exvars.Income(r,t)/exvars.Income(r,t-1),1])^depvars.IncomeElasticityER(r,t,3);
    
    depvars.VTravelDemand(r,t,3) = depvars.VTravelDemand(r,t-1,3)...
        * (exvars.TravelPrice(r,t,3)/exvars.TravelPrice(r,t-1,3))^-(indepvars.PriceElasticity(r,1,3)*2/3)...
        * (exvars.TravelPrice(r,t-1,3)/exvars.TravelPrice(r,t-2,3))^-(indepvars.PriceElasticity(r,1,3)*1/3)...
        * (depvars.EquipmentRate(r,t,3)/depvars.EquipmentRate(r,t-1,3))^indepvars.IncomeElasticity(r,1,3)...
        * (1+indepvars.TimeTrend(r,1,3));
    
    depvars.TravelDemand(r,t,3) = depvars.VTravelDemand(r,t,3) * depvars.EquipmentRate(r,t,3);
    
    % Public Transport (effect by car energy prices!)
        depvars.TravelDemand(r,t,2) = depvars.TravelDemand(r,t-1,2)...
        * (exvars.TravelPrice(r,t,3)/exvars.TravelPrice(r,t-1,3))^(indepvars.PriceElasticity(r,1,2)*2/3)...
        * (exvars.TravelPrice(r,t-1,3)/exvars.TravelPrice(r,t-2,3))^(indepvars.PriceElasticity(r,1,2)*1/3)...
        * (exvars.Income(r,t)/exvars.Income(r,t-1))^indepvars.IncomeElasticity(r,1,2)...
        * (1+indepvars.TimeTrend(r,1,2));
    
    % Air transport
        depvars.TravelDemand(r,t,4) = depvars.TravelDemand(r,t-1,4)...
        * (exvars.TravelPrice(r,t,4)/exvars.TravelPrice(r,t-1,4))^-indepvars.PriceElasticity(r,1,4)...
        * (exvars.Income(r,t)/exvars.Income(r,t-1))^depvars.IncomeElasticity(r,t,4)...
        * (1+indepvars.TimeTrend(r,1,4));
    
end

    depvars.ModeSplit(r,t,modes)         ...
        = depvars.TravelDemand(r,t,modes) /sum(depvars.TravelDemand(r,t,modes),3);
%% Errors
    if t<36 % ensures function is only called with historic data
        depvars.ModeWeight(r,t,modes) = sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,modes)),3);% Ensures not all modes are equal, but that car is not the only important mode
        depvars.ObjectiveFunction(r,t) = sqrt(sum((depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes)).^2 ...
        ./ exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes).^2,3)); % used for oPubTransimization
        depvars.Resid(r,t)             =  depvars.ObjectiveFunction(r,t);
    end

end
%% Write output
    varargout{1} = mean(depvars.ObjectiveFunction(r,time));
if nargout>=2
    varargout{2} = depvars;
end
if nargout >=3
    varargout{3} = indepvars;
end
%% Visualizes minimization process
if cfg.LiveMinimizeVisuals{1} > 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}); % bit clumsy to make the code faster
     GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
end 